function  moments =  simulation_results_output( S,   N_years,  output_file_name, yearly_output  )
% Usage:  simulation_results_output( S, N_years, output_file_name )
% This function generate simulation outputs to the output_file_name, which should be an excel file.
if(nargin<2)
    N_years = 10000;
end

if( isfield( S, 'model_sol' ) )
    model_sol = post_processing(S.model_sol);    
else
    model_sol= post_processing(S);
end
par = model_sol.par;  
shocks = shocks_generation( N_years,  par, 1);  
w_init = 0.2;  lambda_init = 0.1;
if(isfield(par, 'lambda_init'))   % If specifying initial lambda for the simulation
    lambda_init = par.lambda_init;
end
dB_vec = shocks.dB_vec;
dN_vec = shocks.dN_vec;
sim = simulate( w_init, lambda_init,  dB_vec, dN_vec,  N_years, model_sol,  par  );
[ rational_moments,  behavioral_moments ] = simulation_moments( sim, model_sol );

if( par.model_name=="behavioral" )
    moments = behavioral_moments;
else
    moments = rational_moments;
end

underlying_lambda = shocks.lambda_MC;
if( nargin>=4 && yearly_output==true )  % yearly output
    tbl_output = moments.tbl_output(1:12:end,:);
    crisis_vec = zeros( length(tbl_output.T), 1 );
    for( i = 1:length(crisis_vec) )
        crisis_vec(i) = sum(moments.tbl_output.crisis( (12*(i-1)+1):(12*(i-1)+12) ))>0 ; 
    end
    tbl_output = table(  tbl_output.T,  tbl_output.w,   tbl_output.lambda,   tbl_output.lambda_rational,  tbl_output.K,  tbl_output.p,  tbl_output.GDP, tbl_output.productivity  , tbl_output.credit_spread, tbl_output.xK,  tbl_output.dN,    crisis_vec,  tbl_output.rf,  underlying_lambda,   tbl_output.kappap,  moments.RK_monthly,  moments.logK_vec_original , ...
            'VariableNames', {   'T',                   'w',                     'lambda',              'lambda_rational',                     'K',                'p',                    'GDP',       'productivity'               , 'credit_spread',                   'xK',         'dN',    'crisis',  'rf', 'underlying_lambda', 'kappap', 'RK_monthly', 'logK_vec_original'}   );
else
    tbl_output = moments.tbl_output;
    tbl_output = addvars(tbl_output,underlying_lambda);    % add the underlying lambda column
end

if(nargin>=3) % if a file name is specified
    writetable(  tbl_output , output_file_name ,'Sheet', 1, 'WriteVariableNames', true); 
    % % Test
    % T = table([1; 2; 3], {'a'; 'b'; 'c'});
    % writetable(T, '/simulated data/test.xlsx', 'Sheet', 1, 'WriteVariableNames', true);
end

moments.tbl_output = tbl_output;